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Abstract: Since horizon formation in global anti-de Sitter spacetime is dual to thermaliza- 
tion of a conformal field theory on a compact space, whether generic initial data is stable 
or unstable against gravitational collapse is of great interest. We argue that all the known 
stable initial data for massless scalars are dominated by single scalar eigenmodes, specifically 
providing strong numerical evidence consistent with the interpretation that initial data with 
equal energies in two modes collapse on time scales of order the inverse square of the ampli¬ 
tude. We further scan the parameter space for massive scalar field initial data and present 
evidence for a novel class of stable or quasi-stable solutions for massive scalars with energy 
spread through several eigenmodes. 
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1 Introduction 

Through the AdS/CFT correspondence, gravitational physics on (global) d + 1-dimensional 
anti-de Sitter spacetime (AdS^+i) is dual to a conformal field theory (CFT) on M x 5'^“^ 
(the boundary of AdS), with classical gravity valid in the strong coupling regime of the 
CFT. Thermal states of the CFT are dual to black holes in the bulk AdS spacetime,^ so 
thermalization of an initial distribution of energy in the CFT is dual to horizon formation, ie, 
gravitational collapse, in the bulk AdS. It seems natural to expect that even small amounts of 
added energy confined to a compact space generically eventually thermalize, so we are led to a 
perhaps surprising conclusion that generic initial data for matter in AdS should lead to black 
hole formation. On the gravitational side of the correspondence, the reason is that massless 
fields can reach spatial infinity and reflect off the boundary in finite time, so energy cannot 

classically — semi-classically, black holes smaller than the AdS curvature radius are thermodynamically 
unstable to Hawking radiation, so the correct dual to a low-temperature thermal CFT state can be either a 
gravitational radiation gas or a small black hole in equilibrium with such a gas. 
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disperse as it would in asymptotically flat spacetime (similarly, massive fields are confined by 
the effective gravitational potential of AdS). 

From this point of view, there are two natural questions. First, are there classes of 
initial data that are stable against gravitational collapse to a black hole in the bulk AdS 
(especially with measure greater than zero), and what is the dual CFT interpretation of 
this initial data? Second, in the case of gravitational collapse, how long does given initial 
data take to form a horizon — or, in CFT terms, how long does given initial data take to 
thermalize? Strictly speaking, both bulk horizon formation and CFT thermalization take 
infinite boundary time, so we can formulate the question more precisely by asking how long 
it takes for energy to spread through a large number of frequency modes, which typically 
happens just before an approximate horizon forms in the bulk.^ We therefore take bulk 
horizon formation, as we will define approximately in section 2 below, as an indication of 
boundary CFT thermalization as well. This second question is actually easier to answer in 
the unstable case: at low perturbation amplitudes, self-gravitation can only become important 
on time scales of order e“^, where e is the amplitude. In fact, the first question is often phrased 
in terms of stability over times of order e~^. 

In 2011, [1] pioneered numerical work on these questions (see also [2-4]). Specifically, [1] 
presented intriguing numerical evidence that perturbations of AdS are generically unstable 
to black hole formation, at the expected time scale for low amplitudes. The original studies 
of massless scalar field matter have since been followed by studies of complex scalar fields 
[5] and modified theories of gravity [6]. Except when there is a mass gap in the black hole 
spectrum (such as AdSa [7, 8] or AdSs Einstein-Gauss-Bonnet gravity [6]), these numerical 
studies have pointed to instability to horizon formation at arbitrarily small amplitudes for 
generic initial data, which occurs along with a cascade of energy into high frequency modes. 
On the other hand, certain initial data appear stable at low amplitudes, and additional types 
of initial data are the subject of some disagreement (see [9-11]). 

The last several years have seen a simultaneous effort to understand AdS gravitational 
collapse in perturbation theory, particularly with the development of an expansion in terms 
of the scalar eigenmodes on a fixed AdS background [9, 12-14]. Among other symmetries, 
the perturbation theory obeys a scaling law A[t) —)• eA{t/e^) (where A is the coefficient of an 
eigenmode in the expansion), which is also observed in numerical calculations at low amplitude 
[1] and, as suggested above, leads to a universal prediction that horizon formation takes a 
time of order in the perturbative regime regardless of field mass or higher curvature 
terms in the gravitational action [15]. The same scaling symmetry is also apparent in a 
perturbative calculation for a thin-shell scalar held prohle [16]. However, the implications 
of the perturbation theory for stability of generic initial data at arbitrarily small amplitude 
are as yet unclear. Conservation laws of the perturbation theory imply the possibility of 
inverse cascades of energy returning to low-frequency modes, which have been verihed in 

^Again, it is not possible for strict equipartition of a finite energy through an infinite number of modes, so 
we have an approximate notion in mind. 
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numerical simulations (see [9], for example, and our results) and which may be associated 
with stability at low amplitudes. Further, [16, 17] have argued that (quasi-)periodic solutions 
of the perturbation theory should persist as stable solutions of the full theory at arbitrarily 
small amplitudes. On the other hand, [18] have argued for a singularity in the time derivative 
of the phases Bn in generic solutions of the perturbation theory,^ which may be related to 
horizon formation in the full nonlinear theory. As a result, stability, instability, or both may 
be generic in the space of initial data (in the sense of being open sets). 

Given the existing tension, it seems useful to take stock of the characteristics of initial 
data that are agreed to be stable at low amplitudes. From the initial studies, both numeri¬ 
cal and perturbative analysis agree that (nonlinearly modified) scalar eigenmodes are stable 
[1, 20]; these are known as oscillons or, in the complex scalar case, boson stars. ^ In fact, per¬ 
turbations around the oscillon solutions are expected to be stable as well [22], in agreement 
with numerical studies. We note here that all subsequently developed solutions for which 
there is numerical evidence of stability on long time scales are perturbed oscillons in the 
sense that their energy spectrum is dominated by the contributions of a single perturbative 
eigenmode; more precisely, one mode has at least twice the energy of any other individual 
mode. For example, initial data of Gaussian shape and width near the AdS scale i is ap¬ 
parently stable [20, 23, 24] because it is actually dominated by a single scalar eigenmode; 
somewhat wider initial data requires a stronger admixture of other modes and is not quasi¬ 
stable. Likewise, the periodic solutions of [25, 26] are dominated by a single mode as shown 
by their spectra, and, though [9, 27] present a more general ansatz for quasi-periodic solu¬ 
tions to the perturbation theory, the physical solutions presented are all dominated by one 
mode. The highest temperature quasi-periodic solution shown in figure 1 of [27] still has the 
j = 0 mode having approximately twice as much energy as the j = 1 mode. It is also worth 
noting that, in most cases, a single mode dominates the spectrum even more strongly in the 
sense that it has more than 60% of the total energy. Those initial data for which apparently 
stable numerical solutions of the full theory have been presented in the literature are all of 
the latter type. While high temperature quasi-periodic solutions of the perturbation theory 
do exist with considerably less energy in the dominant mode, the behavior of these has not 
to our knowledge been evaluated in the full theory. As the reader may suspect, there are 
also controversial cases: [9, II] hinted at stability of initial data superposing two eigenmodes 
based on a perturbative calculation, which has been criticized by [10]. 

In this paper, we take a numerical approach to the stability of scalar fields in AdS at 
small amplitude, specifically looking at two questions. First, in section 3, we consider massless 
scalars in AdS 4 and provide an independent numerical analysis of the two-mode initial data 
of [9] as well as the superposition of three Gaussians used in [28]. We find that, when the field 
is stable against gravitational collapse over times that grow faster than with decreasing 
amplitude, the energy spectrum is dominated by a single scalar eigenmode as described above. 

^ [19] appeared while this work was in the final stages of preparation and argues the singularity is not present 
in a different gauge, specifically where the time coordinate is boundary time. 

'^In the case of metric perturbations, geons as described in [21] are the equivalent stable eigenmodes. 
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(This condition is similar to the condition of [27] that stable solutions are “close” to a quasi- 
periodic solution of the perturbation theory.) Our calculations also point out difficulties in 
determining the reliability of the perturbative expansion even in the low amplitude regime 
when analytical arguments are not available. 

In section 4, we then turn our attention to the case of scalars with mass /r 7 ^ 0 (corre¬ 
sponding to backgrounds for irrelevant operators in the CFT), which have so far only been 
studied numerically in [28]. (We do not consider tachyonic scalars which are allowed for mass 
squared above the Breitenlohner-Freedman bound [29].) Given that massive fields have dif¬ 
ferent stability properties in asymptotically flat spacetime, this is an important test case for 
stability in AdS. 

As a review, consider the gravitational collapse of massive scalars in asymptotically flat 
spacetime, where an initial configuration can either form a black hole or disperse to infinity 
(ignoring self-interactions that, for example, lead to star birth in astrophysics). Given initial 
data of characteristic width A and <C 1, [30] found that collapse proceeds as for a massless 
scalar, with power-law scaling behavior for the horizon radius near the critical point between 
black hole formation and dispersion [31] (known as a type II transition, see [32]). However, 
for A/i ;§> 1, black holes form only above a finite mass gap (a type I transition). In essence, the 
difference in these two types of transition is whether the scalar potential or gradient energy 
dominates the energy of the initial data. For a massive scalar held in asymptotically AdS 
spacetime with curvature radius £, the different ratios to be considered are A/r, ifi and X/i. 
In principle, the stability properties of the system can change whenever these ratios pass 
through unity. 

We present an overview of the behavior of horizon formation time as a function of these 
ratios in section 4, followed by a more detailed analysis of the behavior at low amplitudes. For 
widths A between the Compton wavelength and AdS scale, we hnd two types of particularly 
interesting behavior, depending on whether the scalar is light or heavy compared to the 
AdS scale. For light fields, we find a discontinuity in the horizon radius at formation as a 
function of the amplitude of initial data, possibly corresponding to a change in efficiency 
of thermalization in the boundary CFT. We analyze this behavior in more detail in section 
5. Of particular interest, for heavy scalars, we find a class of stable (over time scales of 
order e“^) initial data; we present evidence in section 6 that these solutions are qualitatively 
distinct from the single-eigenmode-dominated solutions discussed above in that the energy is 
distributed roughly evenly through several modes. 

We conclude with a discussion of our results and begin here with an overview of our 
methods. 
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2 Methods 


The evolution of a massive scalar field in AdS is governed by the Einstein equations and the 
wave equation, 


Gab + J^Qab = 87 r = 0 , ( 2 . 1 ) 

where the mass of the scalar field (j) is fi. Following [1], we choose a spherically symmetric 
metric ansatz in Schwarzschild-like coordinates 

f2, . 

=-o7—77V ( Ae~‘^^dt^ + A~^dx‘^ + sm^{x/t)dQf '~^) , (2.2) 

cos^(x/t) V / 


where d is the number of spatial dimensions. The areal radius is R{x) = t’tan(x/t'), and we 
henceforth work in units of the AdS scale ^ (ie, i = 1). 

The equations of motion are calculated using a Hamiltonian analysis, as in [33]. The 
evolution equations governing the nonlinear system are 


= Ae"'^n, ^t={Ae-^U 


,-<5t 


n* = 
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(Ae '^tan'^ ^(x)<h)_j 
tan'^“^(x) 
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cos^(x) ’ 


(2.3) 


where H = A is the canonical momentum and = (l>^x is an auxiliary variable. The 

metric functions are determined by 


5^x = — sin(x) cos(x)(n^ + ^h^) 
= (tan(x))‘^~^ 




2 cos 2 (x) 


A = 1 - 2 


sin^(x) M 
{d — 1) tan‘^(x) 


(2.4) 

(2.5) 


where M is the mass function and M{t,x = TTl2)^t = 0. We choose 5{t,x = 0) = 0. 
The linearized system is given by 


d-l 


^,tt — 


sin(x) cos(x) 




d- 


cos^(x) 


:= —L(j) 


( 2 . 6 ) 


(see [34] for the mathematical properties of this system). The eigenfunctions of the operator 
L are given by Jacobi polynomials (see [29] for a review) 


ej{x) = KjCos^^{x)pf^‘^ i’^^/^^^/^)(cos(2x)). 


^ J tXj J j -M. ^ 

The normalization and eigenvalues are given by 


(2.7) 


I 2(2j + A±)j!r(j + A±) 
r(j +V/2)r(j ± Vd2 + 4 ^ 2/2 +T) ’ 


Uj = A± + 2j, A± = - ± -y/d? + 4//2. (2.8) 


In this work we choose A+, corresponding to the normalizable modes according to the inner 
product {f,g):= f{x)g{x)taiii‘^~^{x)dx. 
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As in [ 1 ], we define an “energy per mode” by projecting the field onto eigenmodes. 
Specifically, we define Ilj := (j)j := {4>,ej), and 

:= (tan" {d (tan'^ ^(x)A<l>) —cos ^(x), Cj), so the energy spectrum is 

Ej := - ■ (2-9) 

The total ADM energy is Madm = these projections to study how much 

energy is captured in modes up to some jmax and how energy is transferred between modes 
throughout the evolution. 

As far as we are aware, only initial data (ID) that is either one or more Gaussians in 11 
[1, 3, 5-7, 20, 24, 28], a superposition of two eigenmodes [9, 10], (fake) boson stars [20], or a 
specially constructed time-periodic solution [25] has been studied. We consider four types of 
initial data. We primarily study Gaussian data in H, 

n(t = 0 ,x) = eexp ^- = 0,x) = 0 , ( 2 - 10 ) 

similar to that originally studied in [1]. For a Gaussian pulse of width a, like eq. (2.10), 
we define the wavelength to be A = 2(T. As discussed in the introduction, one goal of this 
work is to explore the interplay among the length scales A, i, and l//i. To this end, we have 
considered several different values of a and 

To explore the universality of our results, we also consider initial data in the form of an 
ingoing pulse 

(^(t = 0 , x) = etan^(x) exp ^-’ ^(t = 0 , x) = 0 (t = 0 , x),x- ( 2 - 11 ) 

This data is more difficult to evolve numerically, specifically as collapse ensues, so we do not 
perform simulations for all the values of a and pL used for the Gaussian initial data in 11. 
Nonetheless, as we present in section 4.3 below, our results appear robust against changes in 
initial data. 

For comparison to existing literature, we also consider two other classes of initial data. 
First, [9] studied two-eigenmode initial data in AdS 4 , specifically, 

4>{t = 0,x) = e(eo(x) -|- Kei(x))/3, (2.12) 


where k is chosen freely. Finally, [28] considered a superposition of three Gaussian wavepackets 
in n 

4(tan(x) — Ri)^\ 

J 

for AdS 4 . We comment on both these classes of initial conditions as they appear in the 
previous literature in addition to studying them ourselves. 


, (),(t = 0, x) = 0, (2.13) 


2 e 


n(t = 0, x) = — tti exp 

TT ^' 


2=1 
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Our numerical methods are similar to those used in [6] and are described in greater detail 
in [24] . We use a fixed grid of 2"^ grid points (with the option to restart at higher resolution), 
mainly with a 4th order Runge-Kutta (RK4) spatial integrator. For some initial data we find 
it necessary to use a fifth-order Dormand-Prince (DoPr) method for the spatial integration. 
The DoPr method is typically only needed when the solution approaches collapse and the 
finer spatial scales need to be resolved more accurately. We terminate the simulation and 
determine that a horizon forms at with k = 7 for /r < 20, while for larger 

masses k = 8. We have found that requiring smaller values of k in some cases requires the 
formation of very narrow troughs in the horizon function A{t, x) which are not resolved by 
the spatial grid, so we have chosen these values to reliably report horizon formation times. A 
discussion of the convergence properties of our code is in appendix A. 

3 Stability of massless scalars in AdS 4 

In this section, we present results of simulations in AdS 4 as a comparison to the previous 
literature. We are particularly concerned with two sets of initial data which have been claimed 
to lead to stable solutions at low amplitude: two-mode initial data and three-Gaussian initial 
data. 

3.1 Two-mode initial data 

First, we consider massless scalars with two-mode initial data as in eq. (2.12) and k = 
3/5, following [9, 11, 27] and [10], which is a point of disagreement between the two sets of 
publications. Specihcally, [9] suggests stability of sufficiently small amplitude two-mode data 
for the length of their simulations, including at e = 0.09, while [10] found horizon formation 
at a long but finite time. We undergo an independent study of this initial data for amplitudes 
e = 0.109,0.09 and 0.08 and find that all three simulations end in collapse. Namely, comparing 
the bottom panel of figure la to figures 3 and 4 of [9] and figure 1 of [10] , shows agreement 
with Bizoh and Rostworowski [10]. Below we provide evidence that the discrepancy is due to 
insufficient resolution in [9, 11]. 

Our simulations also show the same scaling properties first noticed in [1], namely that 
e-^n^ (e^t, X = 0) is an approximate universal function for a given initial field profile as 
long is the simulation is still in the perturbative regime. This scaling symmetry is manifest 
in the improved perturbation theory of [9, 12, 13] (also known as “two-time formalism” or 
TTF), as discussed in the introduction. The close agreement after rescaling, as seen in figure 
la, provides evidence that the numerical solutions are still perturbative until shortly before 
horizon formation. Nonetheless, the perturbative TTF solution diverges from the nonlinear 
numerical solution with the same initial data (including our calculation and those of [9- 
11, 27]). Specifically, in figure 3 of [9], n^(t, x = 0) in the TTF evolution including modes 
up to j = 47 falls below the numerical result by approximately an order of magnitude at 
t ~ 200, while the numerical solution still agrees with [10] and our results. Even the TTF 
solution including modes up to j = 200 [27] diverges strongly from the numerical solution; 
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rescaled time e^t. Bottom panel: = 0). 



panel: evolution of Ej in the lowest six modes. 
Initial data has e = 0.09. 


Figure 1: Results from the evolution of two-mode initial data (2.12) with n = 3/5. 


= 0) is about two orders of magnitude too small for t just below 1100. To be clear, 
it is possible that at smaller e the fully nonlinear evolution could survive to a larger value of 
e^t if n^(t, X = 0) turns over, forming a local maximum; whether all amplitudes collapse after 
two local maxima in n^(t, x = 0) or exhibit a larger number is an open question. 

To gain a better understanding of why the numerical and TTF solutions diverge, we 
perform a spectral decomposition of the energy for e = 0.09. The top panel of figure lb 
shows the total energy in modes up to jmax, while the bottom panel shows the energy in 
the lowest six modes.® We first see a cascade of energy into higher modes, followed by an 
inverse cascade of energy back into the lower modes, in agreement with [9]. There is a 
strong cascade into higher modes just before horizon formation. Figure lb shows that at this 
point about 4% of the energy is in modes greater than j = 47 and almost 2% in j > 256. 
Because higher modes are more sharply peaked around the origin, even a small amount of 
energy in higher modes can cause large differences in the magnitude of n^(x = 0,t) and in 
determining whether or not a collapse occurs. The lesson is that perturbative solutions may 
only be reliable if a large number of eigenmodes are considered, even when the full solution 
is well within the perturbative regime. The correct claim of [11] that the TTF evolution with 
jmax = 4:7 accurately captures the dynamics of the lowest modes is immaterial to the question 
of horizon formation because this is a local question near x = 0 in the small amplitude limit. 
It appears that a small amount of energy transferred to higher modes is sufficient to cause 
gravitational collapse/thermalization. 

We have, however, not yet explained the difference in the non-perturbative numerical 
solutions of [9] and [10] at long times. In principle, one possibility is that these two references 
use different time coordinates (ours matches that of [10]).® For comparison, we perform a 


®We actually plot the normalized energy in the modes, Ej = Ej /Madm ■ 
®Note, though, that [27] uses the same time coordinate we do. 

















(a) Boundary time t vs coordinate time t at res¬ 
olution n = 19. 



(b) Top panel: ADM mass vs time at various 
resolutions. Bottom panel: upper envelope of 
= 0) for the same simulations. 


Figure 2: Coordinate transformation and consistency of our code at several resolutions for 
the e = 0.09 evolution. 


numerical coordinate transformation to the time coordinate t of [9, 11], which is chosen such 
that 5{t,x = 7r/2) = 0 (this is conformal time on the boundary of AdS). The two coordinate 
times are related by 


t = 



exp [—X 


= 7r/2)] dt ; 


(3.1) 


we show the relationship between t, i for late times in figure 2a for the e = 0.09 initial 
data. Because the system is weakly gravitating until collapse, At/At ^ 1.02 for most of the 
simulation; however, as the horizon starts to form, time dilation becomes more significant, and 
the boundary time passes more quickly. We find that collapse is delayed by tjy — tjy ~ 20 for 
the e = 0.09 evolution for resolution and horizon formation threshold given by n = 19. This is 
insufficient to account for the lifetime i ~ 1500 of the simulations of [9].^ Furthermore, time 
dilation only stretches the time axis of hgures 1 and 2b and would not add extra oscillations 
as in figures 3 and 4 of [9]. 

Instead, the difference between [9, 11] and [10] appears entirely due to insufficient reso¬ 
lution in the simulations of [9]. Specifically, the consistency and convergence tests presented 
in [9] were stopped much earlier than the times in the evolution that are in question, the 
longest ending at t ~ 600. As an illustration, the simulation depicted in red in figure 2b has 
grid spacing A ~ 9.6 x 10“®, a higher resolution than used in any of the tests in [11]. The 
rapid decrease in the ADM mass shows that late in the evolution our code also suffers from 
loss of convergence if insufficient resolution is used, and the evolution experiences a similar 

^Technically speaking, a lower threshold for horizon formation results in greater time dilation. However, 
even though we have been unable to locate a discussion of the threshold used in [9], we think this is an unlikely 
source of the error simply because our threshold is quite low. 
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“afterlife” to those of [9, 11]. We have determined that this is due to the largest amplitude 
pulse sharpening and eventually being squeezed between grid points. Figure 2b also shows 
higher resolution simulations where the relative change in ADM mass is AM/M ~ 7 x 10“^ 
for the entirety of the evolution. 

Unfortunately, the convergence tests presented in [11] are insufficient to evaluate whether 
or not the numerical solution is converging to the continuum solution. Typically when per¬ 
forming convergence tests either the number of grid points is increased by a constant amount 
or the resolution is increased by a constant factor. The convergence test presented in figure 1 
of [11] does not do either of these. Rather, what the figure shows is that their code converges 
to the solution with A ~ 1.2 x 10“^ as A asymptotically approaches this value. Meanwhile, 
we have also compared numerical values of n^(t, x = 0 ) with [ 10 ] and obtain the same results 
to several significant figures; an exact match is not expected due to slight differences in the 
algorithms used. The results of our study suggest that the disagreement between [9] and [10] 
would be resolved if [9] used sufficient resolution to be in the convergent regime late in the 
evolution. 

At this point, it is worth assessing the importance of analyzing such apparently well- 
trodden ground (as [11] seemingly wonders). As we stated in the introduction above, one of 
the key questions about gravitational collapse in AdS is what classes of initial data can lead 
to stable evolution at small amplitude. While numerical studies cannot answer this question 
definitively, it is still critical to be careful about whether numerics have actually provided 
evidence of stability or not. In this case, equal-energy two-mode initial data, if stable, would 
present a qualitatively distinct class of stable initial data — all the apparently stable initial 
data for massless scalars is dominated by a single eigenmode. Indeed, some of the authors 
of [9] postulated in [14] that the two-mode initial data is close to a quasi-periodic solution 
(of a type formulated in [9]), but even the quasi-periodic solution closest to their two-mode 
initial data has 70% of its energy in the j = 0 eigenmode. It is also important to realize 
that, while perturbative methods are very powerful, they also suffer from stringent resolution 
requirements (ie, many scalar eigenmodes must be included) because gravitational collapse 
is an essentially local question at low amplitude. In other words, as [27] points out, it is 
the behavior of the very high j tail of the energy spectrum that determines the ultimate 
fate of the system. In summary, while it is impossible at this time to rule out stability of 
equal-energy two-mode initial data as e —)• 0 , all current numerical evidence is consistent with 
instability to horizon formation over times of order precisely due to the contributions of 
high j modes. 

3.2 Multiple-Gaussian initial data 

We also study the multiple Gaussian initial data of eq. (2.13) for massless scalars in AdS 4 , 
which was first considered in [28] with parameters {ai,iTi,i?i} = {1,1/4, 0}, {a 2 ,cr 2 , R 2 } = 
{1/4,1/4, tan( 7 r/ 8 )}, and { 03 , CJ 3 , R 3 } = {1/8,1/4,1}. We have simulated the collapse of this 
initial data and find, as in [28], that it is apparently stable for small amplitudes. To understand 
why this data is stable, we performed a spectral decomposition as a function of time. The top 
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(a) Top panel: sums of the Ej up to jmax- Bot¬ 
tom panel: evolution of the energy in the lowest 
four modes, e = 1. 


(b) Top panel: lP{t,x = 0) (left) and the initial 
profile (right) for e = 1. Bottom panel: tn vs e. 


Figure 3: Results for three-Gaussian initial data (2.13). Below an amplitude of e* = 5.756..., 
we find no horizon formation out to a simulation time of t = 1987. 


panel of figure 3a shows the total energy up to mode jmax, while the bottom panel shows the 
energy in each of the lowest four modes, both for e = 1. We find that approximately 82 per 
cent of the energy is in the lowest mode, essentially independent of time. This suggests that 
this initial data is a perturbation about a single mode, which is known to be stable [1, 20, 25]. 
For the same simulation we plot the upper envelope of n^(t, x = 0) in the top left panel of 
figure 3b. At least for the duration of onr simulation, there is no increase in the Ricci scalar 
at the origin, as is naively expected for stable solutions. 

Our initial data for e = 1 is plotted in the top right panel of hgure 3b and appears to 
match the bottom right panel of figure 9 in [28]. We observe qualitatively similar behavior to 
[28], but our quantitative results differ significantly. For large amplitudes we observe rapid 
horizon formation, but, with decreasing e, we find a small region where the collapse time 
quickly increases and then decreases again. Finally, near e = 5.756 we observe another rapid 
increase in formation time. Yet another small decrease in collapse time is observed followed 
by another increase. It is unclear if this behavior will recur indefinitely. In contrast, [28] 
observes a sudden increase in collapse time at e ~ 1.75 and no earlier increase (see their 
figure 8). 

4 Collapse of massive fields in AdSs 

We now turn to an overview of the behavior of massive scalar fields in AdS. For specificity, 
we work in AdSs, which corresponds to a dual 4D gauge theory. 

Our primary motivation is to consider different values of the dimensionless parameters 
A/U, ifi, and X/i and to determine when the behavior of massive fields diverges from that of 
massless helds. For < 1 (light helds), we consider initial conditions with width satisfying 
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Xfj, < Xji < 1, Xfi < 1 < X/£, and 1 < A/x < X/i. Similarly, for > 1 (heavy fields), we 
consider the X/£ < A/x < 1, X/i < 1 < Xfi, and 1 < X/i < A/x cases. In most cases, we find 
that decreasing initial amplitude leads to monotonically increasing horizon formation times, 
though the behavior may differ in detail from massless scalars. However, as for massless 
scalars, we also find evidence of stability against horizon formation at small amplitudes for 
A ~ when a single eigenmode dominates the spectral distribution. We have also uncovered 
an apparently distinct class of (quasi-)stable initial conditions, which are not dominated by 
any single scalar eigenmode. We discuss those in more detail in section 6. 

We consider multiple classes of initial conditions, demonstrating that the behavior dis¬ 
played is robust. 

4.1 Overview 

For narrow pulses, A less than both l//x and i, we expect the scalar fields to act effectively 
massless, whereas wider pulses may exhibit different behavior. More precisely, since [30, 35] 
found a phase transition at A^ ~ 1, we also may expect a transition in behavior as the initial 
pulse width increases past either the AdS radius or Compton wavelength. However, because 
the scalar field can disperse and then reflect back to the origin (recall that massive fields 
are confined by a gravitational potential in AdS, even though they cannot travel to the AdS 
boundary), the initial distinction in pulse width may be erased after sufficient reflections, 
leading to common behavior at low amplitudes (indeed, this is suggested by the perturbative 
analysis of [15]). The AdS radius also sets the scale of the gravitational potential, so the 
cross-over between light and heavy fields and narrow initial pulses compared to i may also 
change the behavior of the gravitational collapse. Heuristically, we expect the behavior of the 
collapse to be controlled by which type of energy dominates: gradient energy, scalar potential, 
or gravitational potential. 

We consider first light scalars, picking /x = 0.5/i for specificity. To start, we choose 
Gaussian initial data in H as given in equation (2.10) for a variety of widths a. As expected, 
narrow initial pulses (A = 2cj = 0.6i) lead to gravitational collapse behavior similar to that of 
massless fields, as shown in figure 4a. In particular, the horizon formation time forms slightly 
sloped “steps” as a function of pulse amplitude, with each step separated by a jump of Atn 
increasing from approximately 2.8 at large amplitudes to vr at small amplitudes. The initial 
horizon radius xh follows the characteristic arc pattern familiar from [1] with one arc per step 
in tn (though the steps occur rapidly enough at small amplitude that the arcs are difficult 
to resolve). The behavior for pulses wide compared to the AdS scale, shown in figure 4d for 
A = 8i, is similar. The arcs in xh are less distinct at low amplitudes, and more investigation 
would be required to determine the behavior of the horizon radius in this case. Additionally, 
the steps in tn are more steeply sloped, leading to a noticeably larger collapse time even for 
amplitudes that collapse promptly (without reflections from the boundary). At intermediate 
widths, £ < A < 1//X, however, a different structure emerges. Specifically, in the amplitude 
regime where the scalar field undergoes prompt collapse, the initial horizon radius undergoes 
a sudden decrease (not associated with a step in tn). This behavior is apparent in figure 4b 
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(c) xh and tn vs e for a = 0.8£ 



(d) Xh and tn vs e for cr = ii 


Figure 4: Scaling of initial horizon radius xh and formation time tn for mass fj, = 0.5/i 


for a = 0.6£ and in figure 4c for a = O.Sf. We present a further discussion of this phenomenon 
in section 5 below. 

We now turn to the case of heavy scalars; figure 5 shows our results for a held of mass 
H = 5fi. For narrow pulses X < 1/fi, we hnd again, as expected, behavior similar to massless 
helds for the amplitudes we study; results are shown in hgure 5a for a = 0.05£. However, 
the case of heavy scalars and narrow initial data is numerically challenging, likely due to the 
narrow held pulse propagating close the boundary, where the (j)/ cos^(x) term in the scalar 
equation of motion (2.3) becomes subject to large numerical error (since both numerator and 
denominator are small). Wider initial pulses are less numerically challenging. Very wide initial 
pulses (A > £) as in hgure 5d show a lengthened collapse time for prompt collapse surpassing 
tn = 7r/2, which corresponds to the length of time needed for the broadly distributed held 
to collapse toward the center of AdS. The subsequent steps in tn appear more compressed 
as a function of e, leading to a rapid increase in horizon formation time with decreasing 
amplitude. This is apparent to some extent also in hgure 6a for mass ^ = XOjl. We also 
consider intermediate widths Xj< \ < i. AtiT = 0.5£ as in hgure 5c, the step pattern in 
tn has mostly disappeared, replaced by a continuously varying behavior down to a critical 
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amplitude e ss 2.75. (There appear to be steps in tn just above this amplitude, but, since 
the step is /S.tH < vr, it seems that the continuous behavior has simply become steep.) We 
have not been able to form a horizon below this critical amplitude, even while allowing the 
simulations to run to t = 500 (see section 4.2 below for further information). As it turns out, 
we should expect this initial data to be stable; over 91% of the energy of the initial data is 
carried by the j = 0 eigenmode, so this is a single-mode-dominated solution. 

At a somewhat lower widths still in the intermediate range, we find more interesting 
behavior that becomes more apparent at larger masses. For fj, = 5ji, a = 0.3i, we find 
behavior very similar to the massless scalar, at least for solutions with tn ^ 85£. However, 
at the larger mass = and intermediate width a = 0.1^, the horizon formation time tn 
has a sudden narrow jump before decreasing again; see figure 6b. At lower amplitudes, tn 
increases rapidly again (faster than e“^), leading to apparently stable behavior. Note that this 
behavior is distinct from single-mode-dominated solutions as exemplified in figure 5c — in fact, 
as discussed in section 6 below, the energy is distributed democratically throughout several 
eigenmodes. In distinction to the single-mode-dominated oscillon and their quasi-periodic 
generalizations, our results suggest the presence of stable oscillaton solutions for heavy scalar 
fields. For appropriate values of width, then, our initial data can approach the oscillaton 
solutions, leading to an extended collapse time. This effect may also be a manifestation 
of the dynamical mass gap found by [30] in a regime where the mass term dominates over 
gradient energy but the pulse is too narrow to be deformed by the AdS curvature. We revisit 
these issues in more detail in section 6. 

4.2 Low amplitudes 

As mentioned in the introduction, one of the most important questions is the dependence of tn 
on e. At small amplitudes, self-gravitation is suppressed by two powers of e, so gravitational 
collapse should not occur until a time of order e~^ has past. In fact, the scaling symmetry of 
the perturbation theory [9, 12] makes this argument precise; if a solution remains perturbative 
until a time Iq, then a solution with lower amplitude will also remain perturbative until at 
least time to/e^. Since gravitational collapse should require non-perturbative gravitational 
physics, then tn ^ l/^^- While first presented for massless scalars, this argument generalizes 
to massive scalars and Einstein-Gauss-Bonnet gravity [15]. Since it would be noteworthy 
if gravitational collapse happened on a faster time scale (in the perturbative regime), we 
investigate the functional dependence of tn on e in some of our simulations. Similarly, if tn 
grows faster than with decreasing amplitude, that is an indication of quasi-stability (if 
not absolute stability). 

Specihcally, we ht tn = ae^ + b to low amplitude data points, both with and without 
the constraint 6 = 0, and take rough agreement between the values of p in the two fits as an 
indicator that we are in the perturbative regime. We have also checked that the results of 
the fit are robust against removal of some of the low amplitude data points. Results appear 
in table 1 and are generally consistent with tn oc scaling for the initial data shown in the 
table. Due to time constraints, we have not carried out calculations for all our initial data into 
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(b) xh and tn vs e for a = 0.3£ 



(d) Xh and tn vs e for cr = 2i 
Figure 5: Scaling of initial horizon radius xh and formation time tn for mass fi = 5/i 




Figure 6: xh and tu vs e for higher masses 


the perturbative regime; fits to tn from these initial data are characterized by a disagreement 
between the values of p for fits with 6 = 0 and 6 unconstrained. From inspection, it appears 
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/X = 0.5, a = 0.3 

/X = 0.5, (7 = 4 

/X = 5, cr = 2 

^ = 10, (7 = 2 

a 

318 

0.026 

0.077 

0.027 

p 

-2.18 

-2.27 

-2.51 

-2.38 

a 

325 

0.061 

0.446 

0.095 

p 

-2.06 

-2.08 

-1.99 

-2.08 

b 

-7.67 

-13.0 

-25.4 

-20.4 


Table 1: Fits to tn = ae^ + b for indicated parameters. First line constrains 6 = 0. 

that initial data that collapses with tn ^ 100 is typically not yet in the perturbative regime. 

We also noted that scalars with mass fj, = 5ji and initial data width a = 0.5£ lead to 
apparently stable evolution below a critical amplitude e ~ 2.75, and we have been unable 
to induce gravitational collapse below this amplitude. To support this interpretation, we 
have run simulations at amplitudes e = 2.74, 2.67, 2.49 to time t = 500 at a base resolution of 
n = 14. The ADM mass is conserved to 2 parts in 10^^ over this time for all these simulations, 
and we have also run convergence tests (to t = 235, 500, 275 respectively) to verify agreement 
at higher resolution. As we noted above, we expect this initial data to be stable at low 
amplitudes, since it is single-mode dominated, and it remains dominated by the j = 0 mode 
to late times. Specifically, we find that over 80% of the energy remains in the j = 0 mode to 
t = 500£ for all three amplitudes, and the high j spectrum decays exponentially. Interestingly, 
the initial data appears to become stable outside of the perturbative regime — the Ricci scalar 
at the origin does not respect the perturbative scaling symmetry at these three amplitudes, 
and the fraction of the energy held in the j = 0 mode averaged over long times increases with 
decreasing amplitude. 

4.3 Robustness against change of profile shape 

We have also examined the stability behavior of other initial data to determine the robustness 
of results. Specifically, we have considered an ingoing pulse as in eq. (2.11) and two-mode 
initial data (2.12), both in AdSs. 

We considered two values of (/x, a) for comparison of the ingoing pulse to the Gaussian 
initial data in ft, /x = 0.5 with a = 0.3 and /x = 20 with a = 0.1. The initial horizon radii and 
formation times are shown in figure 7 and share the basic characteristics of the corresponding 
mass and width in figures 4a and 6b. In particular, the /x = 20, cr = 0.1 case has the same 
striking sudden increase and reduction in tn as the amplitude decreases. This provides clear 
numerical evidence that our results are robust against change in the shape of initial data. 

Since single mode oscillons and perturbations around them are stable [1, 9], the simplest 
possible data that could produce a horizon is two-mode initial data. We choose the modes to 
have equal energy so they are as far away from a single-mode-dominated solution as possible 
and for comparison with previous work [1, 9] in AdS 4 and [18] in AdSs. With this choice, the 
characteristic width A of the initial data is fixed for a given mass; changing the ratio A/x of 
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(a) = 0.5/£, cr = 0.3^ (b) = 20/£, a = Q.U 


Figure 7: xh and tn vs e for ingoing pulse initial data 


the width to the Compton wavelength requires considering several masses. We consider mass 
values /i = 0, 0.5,1.0 and vTo and estimate the width by htting a Gaussian to the profile (the 
massless case was previously studied in [18]). For ^ = 0.5, A/i < 1; for // = 1, A/i 1; and 
for /r = y/W, A/i > 1. The initial horizon radius xh and horizon formation time tn in these 
cases are qualitatively similar to the results shown in section 4.1 above. 

In particular, we find no evidence that these solutions are close to a stable region. For 
/i = 0,1/2,1 and y/W, we hnd only a direct cascade of energy to higher modes for the duration 
of our simulations, in agreement with the results of [18]. The top panel of hgure 8 shows the 
turbulent energy cascade for the /r = y/W and e = 0.02. While we cannot rule out that inverse 
cascades appear at lower amplitudes, we note that this initial data already appears to be in 
the perturbative regime, as n^(t, x = 0) obeys the perturbative scaling law,® as we see in the 
bottom panel of figure 8. This suggests that AdSs may have smaller stable regions around 
single mode data than AdS 4 , at least for the values of /r we considered. While most of the 
behavior is similar for massive and massless scalars with two-mode initial data, n^(t,x = 0) 
displays one interesting difference: it increases smoothly with small oscillations for massive 
scalars, while for massless scalars it increases in a piecewise linear fashion. 

We also consider the late-time energy spectrum for all the types of initial data we study. 
In collapsing solutions, [24] provides numerical evidence that the spectrum takes the form is 
Ej oc where a = 6/5-1- A[d — 3)/5. We show the late-time spectrum in hgure 9 for several 
masses of two-mode initial data as well as ingoing wave and Gaussian in 11 initial data. These 
data have been evolved to long times {t ~ 1100 for two-mode initial data and t ~ 350 for the 
others). Remarkably the exponent a appears to be independent of both the type of initial 
data and the scalar mass. This suggests that it is the gravitational physics that sets the decay 
rate, rather than the matter content. 


®Note that is not the Ricci scalar at the origin when considering massive scalars, 
verified that the plots are qualitatively identical if one plots the Ricci scalar instead. 


However, we have 
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Figure 8 : Two-mode initial data. Top panel: Energy decomposition for e = 0.02 and 
/r = ^/W. Bottom panel: for several evolutions at different masses and amplitudes. 



Figure 9: Energy spectra late in evolution. Black line shows Ej oc j Initial data are 
two-mode except 11 ID as eq. (2.10) and (j) ID as eq. (2.11). 

5 Initial horizon radius decrease 

In this section we present a more detailed discussion of the rapid decrease in the initial horizon 
radius for length hierarchy ^ < A < l//r, as seen in figures 4b and 4c. As shown in figure 
10a, we have also found similar behavior for a mass /r = 0.1/£ and width a = i. While 
the decrease in xh is not as dramatic in this case, the left-most data point does show a clear 
drop, indicating that this behavior may be common. The astute reader will also notice several 
sudden increases in the initial horizon radius as e decreases in hgures 4b, 4c, and 10a. These 
increases, though again less striking, are familiar from the original work of [1] on massless 
scalar collapse in AdS 4 . 

Both of these phenomena are related to the existence of multiple local minima in the 
metric function A, which we use as a determinant of horizon formation; a typical-looking 
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(a) Discontinuities of xh (b) Horizon function A for /r = 0.5/£, a = 0.8^ 

Figure 10: (a) More evidence of discontinuities in horizon radius for /r = 0.1/£, a = i (b) 

Horizon function A just prior to horizon formation. Red solid curve is e = 1.98, blue dashed 
is e = 1.95. Resolution is n = 14. 


example of this function just prior to tn for ^ = 0.5/£, a = 0.8£, e = 1.98 is shown in 
figure 10b (solid red curve).® Recall that we terminate our calculations and declare horizon 
formation when A decreases below a resolution-dependent threshold; the minimum in A first 
reaches the threshold at When A has multiple local minima (which are associated 

with shells of infalling matter), the initial horizon radius xh is the position of the first local 
minimum to reach the threshold. It is easy to see how decreasing the initial amplitude might 
lead to a jump in xh, then: at some value of e, the innermost shell of matter is dense enough 
to push the inner minimum of A below the threshold, but, at a slightly lower value of e, the 
inner minimum does not contain quite enough mass to decrease past the threshold. Horizon 
formation must wait until another shell of mass approaches the origin, but the jump in mass 
corresponds to a jump in the Schwarzschild radius. It is worth noting that tu increases 
negligibly in this process; the formation of the inner minimum is associated with growth 
in the metric function —6, leading to time dilation that allows the outer mass shells to fall 
inward while very little proper time t at the origin passes. Previous literature [6] has confirmed 
that this behavior remains under increasing resolution, which corresponds to decreasing the 
threshold, but the location of the jump may not be robust as decreasing the threshold can 
cause the jump to shift to higher amplitudes. 

The reverse behavior — a sudden decrease in the horizon radius — is more curious. For a 
fixed resolution (n = 14 in this case), decreasing the amplitude from e = 1.98 to 1.95 changes 
the metric function A from the multi-minimum red solid curve to the single-minimum blue 
dashed curve in figure 10b. Counter-intuitively, the minimum for the lower amplitude is deeper 
than the inner minimum of the higher-amplitude curve. This behavior at fixed amplitude is 
not robust against decreasing the threshold for horizon formation (which, in our analysis, 

® Similar functional forms of A with multiple minima were also found for the case of massless scalar collapse 
in Einstein-Gauss-Bonnet gravity [6]. 
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(a) e = 1.95 (b) e = 1.98 

Figure 11: Energy spectra at horizon formation for ^ = 0.5/£, a = 0.8£ at resolution n = 14. 

occurs with increased resolution); we have verified for e = 1.95 that the initial minimum in 
A (pictured in figure 10b for resolution n = 14) does not pass the threshold for resolution 
n = 16 and that a second local minimum of A appears; the sudden decrease in horizon radius 
should shift to lower amplitudes. However, to be clear, this is a real feature of the solutions 
and not a convergence issue — the n = 14 and higher resolution n = 15,16 solutions agree 
to roundoff, at least until the horizon threshold is reached for n = 14. We have also verified 
convergence of the solution at lower resolution (base resolutions n = 10,12). This is yet 
another way in which gravitational collapse in AdS is sensitive to initial conditions and also 
numerical algorithms. These results also stress that it is crucial to use the same parameters 
when comparing simulations to each other — while most comparisons will not evince this 
type of sensitivity, we have not been able to predict when such sensitivity will appear. In 
vernacular, it is important to compare apples to apples. 

It is at first glance unclear whether we should interpret this phenomenon in the dual 
gauge theory; the initial horizon radius depends on the spatial slicing used to describe the 
gravitational collapse, and nothing physical in the boundary CFT should depend on a gauge 
choice in that way. In addition, the sharp decrease in xh apparently depends sensitively 
on the threshold for horizon formation. However, let us momentarily give in to temptation 
and give our results an interpretation in the boundary theory. Since we take approximate 
horizon formation as corresponding to partial thermalization, we can think of the mass inside 
the initial horizon as corresponding to the fraction of the boundary energy that is approx¬ 
imately thermalized. If we take these results at face value, they tell us that decreasing the 
initial amplitude past a critical value leads to approximate thermalization with less energy 
transferred into higher modes. In figure 11, we display energy spectra for our samples with 
= 0.5/i, a = 0.8.^, and amplitudes e = 1.95 and 1.98 which support this view to some 
degree. Specifically, for e = 1.95, a significantly greater fraction of the energy remains in the 
two lowest modes, indicating that a smaller amount of energy has approximately thermalized. 
Conversely, the energy fraction in higher modes (j > 100, for example) is significantly larger 
for the e = 1.98 evolution. 
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The boundary time of horizon formation provides a possible explanation for this phe¬ 
nomenon. Although both amplitudes lead to similar values of tn as measured in terms of 
proper time at the origin, they have quite different conformal time tn at the boundary due 
to time dilation between the inner minimum of A and the boundary. At resolution n = 14, 
we have tn ~ 0.85 and tn ~ 2.28 for e = 1.95 as compared to tn ~ 0.82 and tu ~ 4.04 for 
e = 1.98. In essence, our measure of thermalization {A passing a fixed threshold) has allowed 
the higher amplitude solution a longer time to thermalize in the boundary theory. Indeed, if 
we increase resolution to n = 16, the e = 1.95 solution develops more minima in A, and the 
energy spectrum 11a evolves to resemble figure 11b more closely. For future work, it would be 
interesting to determine what other simple indicators beside a threshold for A make accurate 
measures of partial thermalization in the boundary theory. 

6 Quasi-stable solutions 

In section 4.1, we noted that massive scalars with intermediate pulse widths 1/^ < A < £ 
exhibit a curious behavior of tn as the amplitude decreases, which we specifically observed 
in figure 6 b for /r = 20 , fi = 0.1 with initial data given by eq. ( 2 . 10 ). At large amplitudes, 
tfj increases in the typical piecewise (roughly) constant fashion. The increase then becomes 
quite rapid, followed by a sudden jump in tn by more than a factor of 3. Very quickly, tn 
decreases again, followed by a steady but rapid increase. This type of behavior has appeared 
in the previous literature for massless [20] and massive [28] scalars. In the massless scalar 
case, the corresponding initial data leads to stable, single-mode-dominated oscillon solutions 
at low amplitudes; in this section, we provide preliminary evidence that the corresponding 
ID for massive scalars may represent a novel class of stable solutions — oscillatons — with 
energy spread democratically through multiple modes. While this behavior also appears for 
ingoing wave initial data as in figure 7b, in this section we will restrict to ID in the form 
( 2 . 10 ) for simplicity. 

At the lowest amplitudes we have studied for = 20, a = 0.1, we hnd that tn rapidly 
increases as e decreases, which provides some support for possible stability at arbitrarily small 
amplitude. These calculations require very high resolution (n > 17) to remain convergent 
until collapse, and we have carried out convergence testing to be conhdent in our calculations. 
In table 2 below, we show fits tn = ae^ + h for the indicated amplitude ranges (which are 
entirely below the large jump in tn noted in section 4); due to computational complexity, 
these fits all have e > 9.62. Except for the smallest amplitude range 9.62 < e < 9.8, we 
find a disagreement in p values for fits constrained to have 6 = 0 or not, and we also find 
that the smaller ranges have increasingly negative values for the power p when 6 = 0. This 
indicates that tn = + 6 does not provide a good fit to tn'-, we have tried also fits of the 

form tn = o(e — e*)^ -|- 6 , but those also fail to provide a good fit (in fact, common fitting 
algorithms do not find sensible fits at all). The lesson of table 2 is not that the horizon 
formation time follows a specific power law in amplitude. What is evident, though, is that tn 
increases much more rapidly than with decreasing amplitude, so this initial data appears 
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e < 11.72 

e < 11.4 

e < 10.6 

e < 9.8 

a 

76000 

1.3 X 10^ 

3.2 X 10^ 

8.2 X 10® 

p 

-2.82 

-3.05 

-3.45 

-4.88 

a 

6.2 X 10^ 

6.1 X 10® 

5.9 X 10® 

7.1 X 10® 

p 

-4.94 

-4.92 

-4.87 

-4.67 

b 

41.4 

38.5 

32.3 

-49.2 


Table 2: Fits to tn = + b for // = 20, a = 0.1 for the indicated amplitude range. First 

line constrains 6 = 0. 




(a) Ej for the lowest six modes. (b) Upper envelope of = 0) 

Figure 12: Evolution with ^ = 20/1', a = 1/10, and e = 9.5. 


to be quasi-stable, that is, stable over longer time scales than expected from generic arguments 
about perturbation theory. Given constraints on computing time, we have not yet ascertained 
whether tn undergoes repeated jumps such as the one seen near e = 11.74 in figure 6b or 
reaches a critical amplitude below which horizons do not (apparently) form. However, it 
seems reasonable to conclude that the low amplitude behavior indicates (quasi-)stability over 
longer than normal time scales at the least. 

We plot the energy in the lowest six modes alongside the upper envelope of n^(t, x = 0) 
in figure 12 for e = 9.5. An inverse energy cascade occurs at t ~ 35 and then again at 
t ~ 130 with a rapid increase in H^ occurring at t ~ 307. There is a pattern of increases 
and decreases in H^ and an approximate recurrence in the spectrum before collapse, both 
consistent with direct and indirect energy cascades. This is similar to the orbits around 
quasi-periodic solutions studied by [27] for massless scalars; however, the spectrum in figure 
12 a is distinct in that the the j = 0,1,2 modes all carry the greatest share of the energy 
at some point in time (compare to the spectrum of equal-energy two-mode ID evolution in 
figure lb). Interestingly, during the final growth in H^ (f ~ 311), ~ 99% of the energy is in 
the lowest 16 modest. 

The presence of inverse cascades tend to coincide with rapid increase in formation time 
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1.736 



Figure 13: Evolution of /x = 20/t', fi = £/10 ID for e = 11.742,11.739,11.736. Top panel: 
upper envelope of TZ{t,x = 0) Bottom panel: Ej for lowest 6 modes for e = 11.739. 

and has been used to argue indirectly for stability in certain solutions (for example, [9]). As 
figure 12b shows, it is possible for collapse to occur even long after an inverse cascade with 
no obvious way of predicting whether or not a rapid cascade of energy to higher modes will 
occur later in the evolution. As we found for two-mode data in section 3.1, the evolution 
may be apparently quasi-periodic for long times and then abruptly change, possibly leading 
to horizon formation. 

Another example of this type of behavior is found in the region near e ~ 11.74, where 
tn jumps from approximately 37.7 to 126 and back to 75.6. We have carried out convergence 
testing on numerical evolution of three amplitudes e = 11.742 (tn 37.73), e = 11.739 
{tn ~ 126.02), and e = 11.736 {tn ~ 78.71) to verify the expected 4th-order convergence. 
The data for the e = 11.739 simulation is presented in appendix A as an example of our 
convergence testing. Furthermore, it is worth noting that we find the long-lived {tn ~ 126.02) 
behavior for more than one amplitude value. As we show in figure 13, evolution at these three 
amplitudes is in remarkable agreement until immediately before horizon formation. While 
that is to some extent to be expected due to the very small changes in amplitude (about 
0.06% over the whole region), it is striking how unpredictable horizon formation seems in a 
comparison of the three amplitudes. In fact, the curvature for the intermediate amplitude, 
which has the largest tn, starts growing but then turns over, allowing a lower amplitude 
evolution to form a horizon hrst. Figure 13 also shows the low j part of the spectrum, which 
evinces a similar pattern of recurrences as the lower amplitude (e = 9.5) evolution shown 
in figure 12a, though the higher-amplitude evolutions have a number of higher-frequency 
oscillations on top of the long-term recurrence pattern. Once again, the spectra for the three 
amplitudes near e = 11.74 are difficult to distinguish except just before one of the evolutions 
forms a horizon. 

To investigate how common this behavior is for massive scalars, we also study extremely 
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massive scalars with [i = 100/£. This data is numerically much more difficult to evolve than 
the less massive fields. Therefore, all simulations discussed here were typically run with of 
214 Qj. more grid points. We then ran simulations at higher resolutions to test if the results 
agreed, as well as did global convergence tests of the longer simulations. We find that the 
solutions that collapse at “late” times (after t ~ 15) require 2^^ or more grid points to be in 
the convergent regime. Running simulations to several hundred time units at this resolution 
is a daunting task, so we present only preliminary results here. See Appendix A for a detailed 
discussion of the convergence of our methods. 

We are interested in data that is not a perturbation about a single mode, so we choose 
a narrow pulse of cr = 1/16; the width of the Gaussian that best fits eo(x) for /r = 100 is 
a K, 7/50 (for // = 20, it is a 3/10). We plot the energy spectrum at t = 0 in the top panel 
of figure 14a. The two lowest modes have nearly equal energy, and an exponential decay in 
the higher modes, so this is not single-mode dominated, at least initially. We also show the 
evolution of the energy in the lowest four modes in the bottom panel of figure 14a. Compared 
to other solutions (such as for ^ = 20 discussed above), we find a strikingly periodic and 
rapid transfer of energy into and out of the lowest mode. Based on a fit pocos(nt + Pi) + P 2 
to we find a frequency approaching ss 2 as the amplitude is decreased. Interestingly, 

we also find that the amplitude of the oscillations decrease for smaller perturbations. 

As we have noted, the stable solutions found so far are perturbations about single mode 
solutions, typically with exponentially decaying spectra at higher modes. In figure 14b, we 
provide some suggestive evidence that this initial data may be stable despite initially having 
nearly equal energies in the j = 0,1 modes (and significant energy in the next two modes). 
We plot the horizon formation time in the bottom panel, and n^(t,x = 0) and TZ{t,x = 0) 
for e = 15 in the top panel. Tracking the rapid increase in formation time has proven 
surprisingly difficult, which means that whether or not this data truly is stable is an open 
question. Nevertheless it is clear that increasing the mass of the scalar results in surprisingly 
different dynamics from what is seen in the massless case. 


use Q to distinguish from the eigenfrequencies uoj. 
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(a) Top panel: Ej{t = 0). Bottom panel: evolu¬ 
tion of Ej in the lowest four modes for e = 15. 
Fit to Eq is po cos(r2t -I- Pi) -I- P 2 - 


(b) Top panel: upper envelope across every seven 
reflections from the boundary for = 0) 

and Tl{t,x = 0). Bottom panel: vs e. 


Figure 14: Evolution for /r = 100, ct = 1/16. The reduced upper envelope is plotted because 
the upper envelopes vary rapidly. 


7 Discussion 

The question of stability of perturbations of AdS against gravitational collapse (either ab¬ 
solutely or on time scales of order for amplitude e) is intriguing both as a problem in 
gravity and through its connection to thermalization in strongly-coupled gauge theories via 
the AdS/CFT correspondence. We have investigated what classes of initial data might lead 
to stable evolution for both massless and massive scalar fields. 

As reviewed in the introduction, it is well-established that perturbations around single 
eigenmode oscillons (or boson stars) are stable, at least on long time scales. We propose, 
therefore, that one criterion for stability is that the energy spectrum is initially dominated by 
a single scalar eigenmode; this is compatible with the conjecture of [27] that stable solutions 
orbit quasi-periodic solutions in configuration space. It is simple to check that most of the 
stable solutions for massless scalars present in extant literature are, in fact, single-mode- 
dominated. We have provided a check that the stable three-Gaussian initial data of [28] 
are single-mode-dominated, as well. In addition, we studied the somewhat controversial 
equal-energy two-mode initial data discussed also in [9-11, 27] and and find that claims of 
stability based on amplitudes studied in the literature so far are premature, though stability 
at smaller amplitudes remains an open question. Specihcally, for the amplitudes of initial 
data studied in the extant literature, according to a fixed approximate measure of horizon 
formation, gravitational collapse occurs at a time proportional to e“^. We conclude that the 
disagreement in the literature over stability of this initial data at small amplitudes is due 
to the use of insufficient resolution in the numerical solutions of [9]. It is also important to 
note that the perturbative TTF solutions require high resolution (in the sense of requiring 
a large number of eigenmodes); even a small fraction of energy in high j modes can trigger 
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gravitational collapse since horizon formation is an essentially local process. Furthermore, 
other than comparison to a numerical solution of the fully nonlinear problem, there is no 
simple diagnostic in the literature for the failure of a perturbative solution, so it is not clear 
that perturbative methods are a reliable gauge of stability even over a fixed time scale. 

Starting in section 4, we present a thorough overview of gravitational collapse for massive 
fields in AdSs. We consider evolution for the six regions of parameter space determined 
by whether each dimensionless ratio ni, \/i is larger or smaller than unity and confirm 
in some examples that our qualitative results for each parameter range are robust against 
changes in the initial data. We also confirm that initial data in most of these parameter 
ranges collapse with tu oc e~^, as expected from perturbation theory. In the process, we 
further discovered two types of novel behavior related to thermalization and stability for 
initial data with intermediate widths. 

First, for scalars with mass less than the AdS scale, we noted, in addition to the more 
familiar sudden increases in the initial horizon radius Xff, also the appearance of sudden 
decreases in Xff with decreasing amplitude. Both phenomena are related to the threshold 
for the metric function A{t, x) used to determine approximate horizon formation and the 
appearance of multiple local minima in A at a fixed time. In some cases, as the amplitude de¬ 
creases, the first-formed minimum of A never reaches the threshold, but subsequently formed 
local minima created by additional infalling matter at larger x can reach the threshold. This 
leads to sudden increases in xh- However, we have shown that the opposite case can oc¬ 
cur; sometimes decreasing the amplitude of the initial data actually allows the first-formed 
local minimum of A to decrease beyond a threshold that it cannot reach at a slightly higher 
amplitude. This puts a spotlight on how approximate horizon formation is defined due to 
the sensitivity of results to the threshold. It also motivates the development of other criteria 
either for approximate horizon formation or thermalization in the dual field theory. 

Our most important result is to provide evidence of a new type of (quasi-)stable field 
profile for massive fields and intermediate width 1/fx < \ < L For /r = 100, a = 1/16, there is 
an apparent critical amplitude, below which we have not been able to find horizon formation. 
For ^ = 20, a = 1/10, we find with decreasing amplitude first the usual increase of tn, then 
a sudden jump from tn ~ 37.7 to tn ~ 126 followed by a decrease to tn ~ 75.6, then a 
rapid increase in tn- While we have not yet ascertained whether there might be a critical 
amplitude, our results do argue for stability over time scales of order e“^. In both cases, the 
energy spectrum is not dominated by a single mode, unlike for other known stable solutions, 
but rather is initially spread democratically through several modes, possibly indicating the 
existence of a quasi-stable multi-mode oscillaton solution. Following the hypothesis [27] that 
stable evolutions are orbits of quasi-periodic solutions, it would be interesting to determine 
if quasi-periodic solutions of the perturbation theory for a massive scalar allow energy to 
be distributed more evenly between modes. An alternate possibility is that the “island of 
stability” for perturbations around single-mode oscillons or quasi-periodic solutions is much 
wider for heavy scalars in an AdS manifestation of the mass gap found in asymptotically flat 
spacetime [30]. Another important point is that our fj, = 20, a = 1/10 solutions in the region 
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when tfi jumps from 37.7 —)• 126 —75.6 show a significant agreement up until the time 
of collapse, even though the final horizon formation time differs greatly. There is no clear 
indication in advance whether a given amplitude will collapse sooner or later.We leave a 
more careful study of these potentially novel stable solutions for future work. 

In summary, there are two major lessons from our work. First is that, while there is 
by now a well-developed leading-order perturbation theory for scalar gravitational collapse 
in AdS, comparison to numerical solutions of the full nonlinear theory shows that it is very 
difficult to predict when the perturbative theory will break down. We have given several 
explicit examples of this difficulty; as a result, it is difficult to know how long a perturbative 
solution remains reliable (for example in the sense of [17]). Second, all the known stable or 
quasi-stable evolutions of massless scalars in AdS are dominated by a single eigenmode of 
the scalar on a fixed AdS background, whereas quasi-stable evolutions of the massive scalar 
can have energy democratically spread through several eigenmodes. These massive scalar 
solutions call for future investigation. 
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A Convergence tests and code validation 

As mentioned at the end of section 2, in addition to the standard RK4 spatial integration, 
we also use a fifth order DoPr method. The DoPr reduces the resolution needed to bring 
massive solutions to collapse. When using RK4, we often need a spatial resolution of 2^® 
grid points near the end of long simulations (or those near critical points) to prevent crashes. 
On the other hand, with DoPr, it often suffices to use a resolution of 2^® grid points, but 
DoPr is significantly slower than RK4 at the same resolution. We have also tried different 
spatial integrators. To validate our code, we have checked that several different methods give 
consistent results over a long simulation of equal-energy two-mode ID with e = 0.109, as 
discussed in section 3.1. As shown in figure 15, these methods show a remarkable agreement 
for lP{t,x = 0). Furthermore, they agree on tn with a relative error of 4 x 10“® and on xh 
to within one grid point. 

re-emphasize, these results have been validated by convergence testing and show the expected 4th order 
convergence up to collapse. 
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Figure 15: = 0) as in figure la with e = 0.109 for several numerical methods 


There are two quantities we monitor for consistency of the solutions. The system (2.3) has 
a constraint, C := = 0, and the ADM mass should be conserved by the time evolution. 

Since explicit convergence tests for each simulation are too demanding, we typically monitor 
the ADM mass and C for consistency. For small and zero this is usually sufficient because 
we find that the mass rapidly decreases as convergence is lost (see figure 2b). This happens 
because the pulses slip between the mesh if insufficient resolution is used. However, as /r is 
increased the ADM mass is conserved well and C will remain small even when convergence 
is lost. A careful study of the evolution when this occurs shows that it is the formation of 
several small minima in A that are under resolved. These minima each reach a different 
value and the narrower, farther out (radially) ones are ultimately those that trigger horizon 
formation. We find that in many cases at least 2^^ grid points are required to resolve these 
minima, making the computations incredibly expensive. Nevertheless for particularly long or 
interesting simulations, we have run explicit convergence tests to provide credibility for our 
results. 

We perform an array of convergence tests to determine the reliability of our results, which 
are presented here for jjL = 20, cr = 0.1, and amplitude e = 11.74. This solution has a large 
tfj compared to slightly larger or smaller amplitudes, so it is an interesting test case for 
convergence and consistency of the code. 

Figure 16a demonstrates that the scalar field, mass function, and metric functions A and 
5 are fourth-order convergent (in norm) through multiple reflections off the outer bound¬ 
ary, and figure 16b shows that the constraint residual C decreases with increasing resolution. 
Figure 16c shows that the ADM mass is well-conserved for the duration of our simulation. We 
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(a) Global order of convergence 
of (j), M, A, S 


(b) Constraint violations C (c) ADM mass and Richardson 

extrapolated mass for coarsest 
resolution 



lime time lime 


(d) Convergence of ADM mass 


(e) Order of convergence of (f) Order of convergence of con- 
lP{x = 0,t) straint violation C. 


Figure 16: Several measures of the accuracy and convergence of our numerical methods. The 
base resolution is 2^^ + 1 with refinement by factors of two. The test was performed using 
the DoPr spatial integrator. The initial data used eq. (2.10) with /x = 20, cr = 0.1, e = 11.74. 


plot the ADM mass at all three resolutions as well as the Richardson extrapolated value. 
The order of convergence of the ADM mass is plotted in figure 16d and exhibits approximately 
sixth order convergence. Because of the variable step size method used with the DoPr inte¬ 
gration, the ADM mass has a higher order of convergence than expected solely from the time 
evolution. Since n^(t,x = 0) has been used to detect the onset of the turbulent instability. 
We show convergence of this quantity in figure 16e. Finally, we plot the convergence of the 
constraint residual in figure 16f, which demonstrates that the during the evolution we stay 
on the constraint sub-manifold. In post processing we used a second order method, unlike 
in the actual time evolution where we use a fourth order method. The results of our testing 
demonstrate that our code is convergent and consistent even for long simulations with several 
reflections off the outer boundary, as well as for large scalar field mass values. 

^^See [36] for a discussion of this technique as a measure of consistency with adaptive mesh refinement. 
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